!polyfits
! E02ACF
! E02ADF
! E02ADZ
! E02AEF
! E02AFF
! E02AGF

      SUBROUTINE E02ACF(X,Y,N,AA,M1,REF)
C     MARK 1 RELEASE.  NAG COPYRIGHT 1971
C     MARK 4.5 REVISED
C     MARK 5C REVISED
C     MARK 9B REVISED. IER-361 (JAN 1982)
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C     MARK 14C REVISED. IER-877 (NOV 1990).
C     CALCULATES A MINIMAX POLYNOMIAL FIT TO A SET OF DATA POINTS
C     AS A
C     SERIES OF CHEBYSHEV POLYNOMIALS.
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='E02ACF')
C     .. Scalar Arguments ..
      DOUBLE PRECISION  REF
      INTEGER           M1, N
C     .. Array Arguments ..
      DOUBLE PRECISION  AA(M1), X(N), Y(N)
C     .. Local Scalars ..
      DOUBLE PRECISION  ABSHI, AI, AI1, D, DENOM, H, HI, HIMAX, HMAX,
     *                  ONE, P1, P5, PREVH, RHI, RHI1, XI, XJ, XNEXTH,
     *                  ZERO
      INTEGER           I, I1, I2, IFAIL, IJ1, IMAX, IRI, IRJ, J, J1, K,
     *                  M, M2
C     .. Local Arrays ..
      DOUBLE PRECISION  A(100), RH(100), RX(100)
      INTEGER           IR(100)
      CHARACTER*1       P01REC(1)
C     .. External Functions ..
      INTEGER           P01ABF
      EXTERNAL          P01ABF
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, DBLE, INT
C     .. Data statements ..
      DATA              ZERO/0.0D0/, ONE/1.0D0/, P5/0.5D0/, P1/0.1D0/
C     .. Executable Statements ..
C     ENFORCED HARD FAIL FOR OUT-OF-BOUNDS PARAMETERS
      IF (M1.GE.N .OR. M1.GE.100) IFAIL = P01ABF(0,1,SRNAME,0,P01REC)
      DO 20 I = 2, N
         IF (X(I).LE.X(I-1)) IFAIL = P01ABF(0,2,SRNAME,0,P01REC)
   20 CONTINUE
      M2 = M1 + 1
      M = M1 - 1
      PREVH = -ONE
      IR(1) = 1
      IR(M2) = N
      D = DBLE(N-1)/DBLE(M1)
      H = D
      IF (M.EQ.0) GO TO 60
      DO 40 I = 2, M1
         IR(I) = INT(H+P5) + 1
         H = H + D
   40 CONTINUE
   60 H = -ONE
      DO 80 I = 1, M2
         IRI = IR(I)
         RX(I) = X(IRI)
         A(I) = Y(IRI)
         RH(I) = -H
         H = -H
   80 CONTINUE
      DO 120 J = 1, M1
         I1 = M2
         AI1 = A(I1)
         RHI1 = RH(I1)
         I = M2
  100    I = I - 1
         IJ1 = I - J + 1
         DENOM = RX(I1) - RX(IJ1)
         AI = A(I)
         RHI = RH(I)
         A(I1) = (AI1-AI)/DENOM
         RH(I1) = (RHI1-RHI)/DENOM
         I1 = I
         AI1 = AI
         RHI1 = RHI
         IF (I-J) 120, 120, 100
  120 CONTINUE
      H = -A(M2)/RH(M2)
      DO 140 I = 1, M2
         A(I) = A(I) + RH(I)*H
  140 CONTINUE
      IF (M.EQ.0) GO TO 200
      J = M1
  160 J = J - 1
      XJ = RX(J)
      I = J
      AI = A(I)
      J = J + 1
      DO 180 I1 = J, M1
         AI1 = A(I1)
         A(I) = AI - XJ*AI1
         AI = AI1
         I = I1
  180 CONTINUE
      J = J - 1
      IF (J-1) 200, 200, 160
  200 CONTINUE
      HMAX = ABS(H)
      IF (HMAX.GT.PREVH) GO TO 220
      A(M2) = -HMAX
      GO TO 480
  220 A(M2) = HMAX
      PREVH = HMAX
      IMAX = IR(1)
      HIMAX = H
      J = 1
      IRJ = IR(J)
      DO 300 I = 1, N
         IF (I.EQ.IRJ) GO TO 280
         XI = X(I)
         HI = ZERO
         K = M2
  240    K = K - 1
         HI = HI*XI + A(K)
         IF (K-1) 260, 260, 240
  260    HI = HI - Y(I)
         ABSHI = ABS(HI)
         IF (ABSHI.LE.HMAX) GO TO 300
         HMAX = ABSHI
         HIMAX = HI
         IMAX = I
         GO TO 300
  280    IF (J.GE.M2) GO TO 300
         J = J + 1
         IRJ = IR(J)
  300 CONTINUE
      IF (IMAX.EQ.IR(1)) GO TO 480
      DO 320 I = 1, M2
         IF (IMAX.LT.IR(I)) GO TO 340
  320 CONTINUE
      I = M2
  340 I2 = INT(DBLE(I)*P5)
      I2 = I - 2*I2
      XNEXTH = H
      IF (I2.EQ.0) XNEXTH = -H
      IF (HIMAX*XNEXTH.LT.0.0D0) GO TO 360
      IR(I) = IMAX
      GO TO 60
  360 IF (IMAX.GE.IR(1)) GO TO 420
      J1 = M2
      J = M2
  380 J = J - 1
      IR(J1) = IR(J)
      J1 = J
      IF (J-1) 400, 400, 380
  400 IR(1) = IMAX
      GO TO 60
  420 IF (IMAX.LE.IR(M2)) GO TO 460
      J = 1
      DO 440 J1 = 2, M2
         IR(J) = IR(J1)
         J = J1
  440 CONTINUE
      IR(M2) = IMAX
      GO TO 60
  460 IR(I-1) = IMAX
      GO TO 60
  480 CONTINUE
      DO 500 I = 1, M1
         AA(I) = A(I)
  500 CONTINUE
      REF = A(M2)
      RETURN
      END

      SUBROUTINE E02ADF(M,KPLUS1,NROWS,X,Y,W,WORK1,WORK2,A,S,IFAIL)
C
C     NAG LIBRARY SUBROUTINE  E02ADF
C
C     E02ADF  COMPUTES WEIGHTED LEAST-SQUARES POLYNOMIAL
C     APPROXIMATIONS TO AN ARBITRARY SET OF DATA POINTS.
C
C     FORSYTHE-CLENSHAW METHOD WITH MODIFICATIONS DUE TO
C     REINSCH AND GENTLEMAN.
C
C     USES NAG LIBRARY ROUTINE  P01AAF.
C     USES BASIC EXTERNAL FUNCTION  SQRT.
C
C     STARTED - 1973.
C     COMPLETED - 1976.
C     AUTHOR - MGC AND JGH.
C
C     WORK1  AND  WORK2  ARE WORKSPACE AREAS.
C     WORK1(1, R)  CONTAINS THE VALUE OF THE  R TH  WEIGHTED
C     RESIDUAL FOR THE CURRENT DEGREE  I.
C     WORK1(2, R)  CONTAINS THE VALUE OF  X(R)  TRANSFORMED
C     TO THE RANGE  -1  TO  +1.
C     WORK1(3, R)  CONTAINS THE WEIGHTED VALUE OF THE CURRENT
C     ORTHOGONAL POLYNOMIAL (OF DEGREE  I)  AT THE  R TH
C     DATA POINT.
C     WORK2(1, J)  CONTAINS THE COEFFICIENT OF THE CHEBYSHEV
C     POLYNOMIAL OF DEGREE  J - 1  IN THE CHEBYSHEV-SERIES
C     REPRESENTATION OF THE CURRENT ORTHOGONAL POLYNOMIAL
C     (OF DEGREE  I).
C     WORK2(2, J)  CONTAINS THE COEFFICIENT OF THE CHEBYSHEV
C     POLYNOMIAL OF DEGREE  J - 1  IN THE CHEBYSHEV-SERIES
C     REPRESENTATION OF THE PREVIOUS ORTHOGONAL POLYNOMIAL
C     (OF DEGREE  I - 1).
C
C     NAG COPYRIGHT 1975
C     MARK 5 RELEASE
C     MARK 6 REVISED  IER-84
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     CHECK THAT THE VALUES OF  M  AND  KPLUS1  ARE REASONABLE
C
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='E02ADF')
C     .. Scalar Arguments ..
      INTEGER           IFAIL, KPLUS1, M, NROWS
C     .. Array Arguments ..
      DOUBLE PRECISION  A(NROWS,KPLUS1), S(KPLUS1), W(M), WORK1(3,M),
     *                  WORK2(2,KPLUS1), X(M), Y(M)
C     .. Local Scalars ..
      DOUBLE PRECISION  ALPIP1, BETAI, BJ, BJP1, BJP2, CI, D, DF, DI,
     *                  DIM1, DJ, EPSR, FACTOR, PIJ, SIGMAI, WRPR,
     *                  WRPRSQ, X1, XCAPR, XM
      INTEGER           I, IERROR, IPLUS1, IPLUS2, J, JPLUS1, JPLUS2,
     *                  JREV, K, MDIST, R
C     .. Local Arrays ..
      CHARACTER*1       P01REC(1)
C     .. External Functions ..
      INTEGER           P01ABF
      EXTERNAL          P01ABF
C     .. Intrinsic Functions ..
      INTRINSIC         SQRT
C     .. Executable Statements ..
      IERROR = 4
      IF (KPLUS1.LT.1 .OR. M.LT.KPLUS1) GO TO 380
      K = KPLUS1 - 1
C
C     TEST THE VALIDITY OF THE DATA.
C
C     CHECK THAT THE WEIGHTS ARE STRICTLY POSITIVE.
C
      IERROR = 1
      DO 20 R = 1, M
         IF (W(R).LE.0.0D0) GO TO 380
   20 CONTINUE
C
C     CHECK THAT THE VALUES OF  X(R)  ARE NON-DECREASING AND
C     DETERMINE
C     THE NUMBER  (MDIST)  OF DISTINCT VALUES OF  X(R).
C
      IERROR = 2
      MDIST = 1
      DO 40 R = 2, M
         IF (X(R).LT.X(R-1)) GO TO 380
         IF (X(R).EQ.X(R-1)) GO TO 40
         MDIST = MDIST + 1
   40 CONTINUE
C
C     IF THE  X(R)  ALL HAVE THE SAME VALUE, I.E.  MDIST = 1,
C     THE NORMALIZATION OF THE INDEPENDENT VARIABLE IS NOT
C     POSSIBLE.
C
      IERROR = 3
      IF (MDIST.EQ.1) GO TO 380
C
C     IF THE NUMBER OF DISTINCT VALUES OF  X(R)  FAILS TO EXCEED
C     THE MAXIMUM DEGREE  K  THERE IS NO UNIQUE POLYNOMIAL
C     APPROXIMATION OF THAT DEGREE.
C
      IERROR = 4
      IF (MDIST.LE.K) GO TO 380
C
C     CHECK THAT  NROWS  HAS BEEN SET SUFFICIENTLY LARGE.
C
      IERROR = 5
      IF (NROWS.LT.KPLUS1) GO TO 380
      IERROR = 0
C
      X1 = X(1)
      XM = X(M)
      D = XM - X1
C
C     THE INITIAL VALUES  WORK1(1, R)  (R = 1, 2, ..., M)  OF THE
C     WEIGHTED RESIDUALS AND THE VALUES  WORK1(2, R)  (R = 1, 2,
C     ..., M)
C     OF THE NORMALIZED INDEPENDENT VARIABLE ARE COMPUTED.  NOTE
C     THAT
C     WORK1(2, R)  IS COMPUTED FROM THE EXPRESSION BELOW RATHER
C     THAN THE
C     MORE NATURAL FORM
C
C     (2.0*X(R) - X1 - XM)/D,
C
C     SINCE THE FORMER GUARANTEES THE COMPUTED VALUE TO DIFFER FROM
C     THE TRUE VALUE BY AT MOST  4.0*MACHINE ACCURACY,  WHEREAS THE
C     LATTER HAS NO SUCH GUARANTEE.
C
      DO 60 R = 1, M
         WORK1(1,R) = W(R)*Y(R)
         WORK1(2,R) = ((X(R)-X1)-(XM-X(R)))/D
   60 CONTINUE
      I = 1
      BETAI = 0.0D0
      DO 360 IPLUS1 = 1, KPLUS1
C
C        SET STARTING VALUES FOR DEGREE  I.
C
         IPLUS2 = IPLUS1 + 1
         IF (IPLUS1.EQ.KPLUS1) GO TO 100
         DO 80 JPLUS1 = IPLUS2, KPLUS1
            A(IPLUS1,JPLUS1) = 0.0D0
   80    CONTINUE
         WORK2(1,IPLUS2) = 0.0D0
         WORK2(2,IPLUS2) = 0.0D0
  100    ALPIP1 = 0.0D0
         CI = 0.0D0
         DI = 0.0D0
         A(I,IPLUS1) = 0.0D0
         WORK2(1,IPLUS1) = 1.0D0
         IF (KPLUS1.GT.1) WORK2(2,1) = WORK2(1,2)
         DO 260 R = 1, M
            XCAPR = WORK1(2,R)
C
C           THE WEIGHTED VALUE  WORK1(3, R)  OF THE ORTHOGONAL
C           POLYNOMIAL OF
C           DEGREE  I  AT  X = X(R)  IS COMPUTED BY RECURRENCE FROM ITS
C           CHEBYSHEV-SERIES REPRESENTATION.
C
            IF (IPLUS1.GT.1) GO TO 120
            WRPR = W(R)*0.5D0*WORK2(1,1)
            WORK1(3,R) = WRPR
            GO TO 240
  120       J = IPLUS2
            IF (XCAPR.GT.0.5D0) GO TO 200
            IF (XCAPR.GE.-0.5D0) GO TO 160
C
C           GENTLEMAN*S MODIFIED RECURRENCE.
C
            FACTOR = 2.0D0*(1.0D0+XCAPR)
            DJ = 0.0D0
            BJ = 0.0D0
            DO 140 JREV = 1, I
               J = J - 1
               DJ = WORK2(1,J) - DJ + FACTOR*BJ
               BJ = DJ - BJ
  140       CONTINUE
            WRPR = W(R)*(0.5D0*WORK2(1,1)-DJ+0.5D0*FACTOR*BJ)
            WORK1(3,R) = WRPR
            GO TO 240
C
C           CLENSHAW*S ORIGINAL RECURRENCE.
C
  160       FACTOR = 2.0D0*XCAPR
            BJP1 = 0.0D0
            BJ = 0.0D0
            DO 180 JREV = 1, I
               J = J - 1
               BJP2 = BJP1
               BJP1 = BJ
               BJ = WORK2(1,J) - BJP2 + FACTOR*BJP1
  180       CONTINUE
            WRPR = W(R)*(0.5D0*WORK2(1,1)-BJP1+0.5D0*FACTOR*BJ)
            WORK1(3,R) = WRPR
            GO TO 240
C
C           REINSCH*S MODIFIED RECURRENCE.
C
  200       FACTOR = 2.0D0*(1.0D0-XCAPR)
            DJ = 0.0D0
            BJ = 0.0D0
            DO 220 JREV = 1, I
               J = J - 1
               DJ = WORK2(1,J) + DJ - FACTOR*BJ
               BJ = BJ + DJ
  220       CONTINUE
            WRPR = W(R)*(0.5D0*WORK2(1,1)+DJ-0.5D0*FACTOR*BJ)
            WORK1(3,R) = WRPR
C
C           THE COEFFICIENT  CI  OF THE  I TH  ORTHOGONAL POLYNOMIAL AND
C           THE
C           COEFFICIENTS  ALPIP1  AND  BETA I  IN THE
C           THREE-TERM RECURRENCE RELATION FOR THE ORTHOGONAL
C           POLYNOMIALS ARE COMPUTED.
C
  240       WRPRSQ = WRPR**2
            DI = DI + WRPRSQ
            CI = CI + WRPR*WORK1(1,R)
            ALPIP1 = ALPIP1 + WRPRSQ*XCAPR
  260    CONTINUE
         CI = CI/DI
         IF (IPLUS1.NE.1) BETAI = DI/DIM1
         ALPIP1 = 2.0D0*ALPIP1/DI
C
C        THE WEIGHTED RESIDUALS  WORK1(1, R)  (R = 1, 2, ..., M)  FOR
C        DEGREE  I  ARE COMPUTED, TOGETHER WITH THEIR SUM OF SQUARES,
C        SIGMAI.
C
         SIGMAI = 0.0D0
         DO 280 R = 1, M
            EPSR = WORK1(1,R) - CI*WORK1(3,R)
            WORK1(1,R) = EPSR
            SIGMAI = SIGMAI + EPSR**2
  280    CONTINUE
C
C        THE ROOT MEAN SQUARE RESIDUAL  S(I + 1)  FOR DEGREE  I  IS
C        THEORETICALLY UNDEFINED IF  M = I + 1  (THE CONDITION FOR THE
C        POLYNOMIAL TO PASS EXACTLY THROUGH THE DATA POINTS).  SHOULD
C        THIS
C        CASE ARISE THE R.M.S. RESIDUAL IS SET TO ZERO.
C
         IF (IPLUS1.GE.M) GO TO 300
         DF = M - IPLUS1
         S(IPLUS1) = SQRT(SIGMAI/DF)
         GO TO 320
  300    S(IPLUS1) = 0.0D0
C
C        THE CHEBYSHEV COEFFICIENTS  A(I + 1, 1), A(I + 1, 2), ...,
C        A(I + 1, I + 1)  TOGETHER WITH THE COEFFICIENTS
C        WORK2(1, 1), WORK2(1, 2), ..., WORK2(1, I + 1),   IN THE
C        CHEBYSHEV-SERIES REPRESENTATION OF THE  I TH  ORTHOGONAL
C        POLYNOMIAL ARE COMPUTED.
C
  320    DO 340 JPLUS1 = 1, IPLUS1
            JPLUS2 = JPLUS1 + 1
            PIJ = WORK2(1,JPLUS1)
            A(IPLUS1,JPLUS1) = A(I,JPLUS1) + CI*PIJ
            IF (JPLUS2.GT.KPLUS1) GO TO 380
            WORK2(1,JPLUS1) = WORK2(1,JPLUS2) + WORK2(2,JPLUS1) -
     *                        ALPIP1*PIJ - BETAI*WORK2(2,JPLUS2)
            WORK2(2,JPLUS2) = PIJ
  340    CONTINUE
         DIM1 = DI
         I = IPLUS1
  360 CONTINUE
  380 IF (IERROR) 400, 420, 400
  400 IFAIL = P01ABF(IFAIL,IERROR,SRNAME,0,P01REC)
      RETURN
  420 IFAIL = 0
      RETURN
      END

      SUBROUTINE E02ADZ(MFIRST,MLAST,MTOT,KPLUS1,NROWS,KALL,NDV,X,Y,W,
     *                  XMIN,XMAX,INUP1,NU,WORK1,WORK2,A,S,SERR,EPS,
     *                  IFAIL)
C     MARK 7 RELEASE. NAG COPYRIGHT 1978.
C     MARK 8 REVISED. IER-228 (APR 1980).
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     E02ADZ  COMPUTES WEIGHTED LEAST-SQUARES POLYNOMIAL
C     APPROXIMATIONS TO AN ARBITRARY SET OF DATA POINTS,
C     WITH, IF REQUIRED, SEVERAL SETS OF VALUES OF THE
C     DEPENDENT VARIABLE.
C
C     FORSYTHE-CLENSHAW METHOD WITH MODIFICATIONS DUE TO
C     REINSCH AND GENTLEMAN.
C
C     STARTED - 1973.
C     COMPLETED - 1978.
C     AUTHORS - MGC AND GTA.
C
C     WORK1  AND  WORK2  ARE WORKSPACE AREAS.
C     WORK1(1, R)  CONTAINS THE VALUE OF  X(R)  TRANSFORMED
C     TO THE RANGE  -1  TO  +1.
C     WORK1(2, R)  CONTAINS THE WEIGHTED VALUE OF THE CURRENT
C     ORTHOGONAL POLYNOMIAL (OF DEGREE  I)  AT THE  R TH
C     DATA POINT.
C     WORK2(1, J)  CONTAINS THE COEFFICIENT OF THE CHEBYSHEV
C     POLYNOMIAL OF DEGREE  J - 1  IN THE CHEBYSHEV-SERIES
C     REPRESENTATION OF THE CURRENT ORTHOGONAL POLYNOMIAL
C     (OF DEGREE  I).
C     WORK2(2, J)  CONTAINS THE COEFFICIENT OF THE CHEBYSHEV
C     POLYNOMIAL OF DEGREE  J - 1  IN THE CHEBYSHEV-SERIES
C     REPRESENTATION OF THE PREVIOUS ORTHOGONAL POLYNOMIAL
C     (OF DEGREE  I - 1).
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  XMAX, XMIN
      INTEGER           IFAIL, INUP1, KALL, KPLUS1, MFIRST, MLAST, MTOT,
     *                  NDV, NROWS
C     .. Array Arguments ..
      DOUBLE PRECISION  A(NDV,NROWS,KPLUS1), EPS(NDV,MLAST), NU(INUP1),
     *                  S(NDV,KPLUS1), SERR(KPLUS1), W(MLAST),
     *                  WORK1(2,MTOT), WORK2(2,KPLUS1), X(MLAST),
     *                  Y(NDV,MLAST)
C     .. Local Scalars ..
      DOUBLE PRECISION  ALPIP1, BETAI, BJ, BJP1, BJP2, CIL, D, DF, DI,
     *                  DIM1, DJ, EPSLR, FACTOR, PIJ, SIGMAI, WR, WRPR,
     *                  WRPRSQ, X1, XCAPR, XM
      INTEGER           I, IERROR, II, IM1, INU, IPLUS1, IPLUS2, J,
     *                  JPLUS1, JPLUS2, JREV, K, L, M, MDIST, MR, R
      LOGICAL           WNZ
C     .. Local Arrays ..
      DOUBLE PRECISION  CI(10)
C     .. Intrinsic Functions ..
      INTRINSIC         SQRT
C     .. Executable Statements ..
      K = KPLUS1 - 1
      INU = INUP1 - 1
C
C     TEST THE VALIDITY OF THE DATA.
C
C     CHECK INPUT PARAMETERS.
C
      M = MLAST - MFIRST + 1
      I = KPLUS1 - INU
      IERROR = 5
      IF (MFIRST.LT.1 .OR. INUP1.LT.1 .OR. KPLUS1.LT.INUP1 .OR. M.LT.
     *    I .OR. NDV.LT.1 .OR. (KALL.NE.1 .AND. KALL.NE.0)) GO TO 600
C
C     CHECK THAT THE VALUES OF X(R) ARE NON-DECREASING AND
C     DETERMINE THE NUMBER (MDIST) OF DISTINCT VALUES OF X(R)
C     WITH NON-ZERO WEIGHT
C
      IERROR = 2
      MDIST = 1
      IF (W(MFIRST).EQ.0.0D+0) MDIST = 0
      L = MFIRST + 1
      IF (L.GT.MLAST) GO TO 40
      WNZ = W(MFIRST) .NE. 0.0D+0
      DO 20 R = L, MLAST
         IF (X(R).LT.X(R-1)) GO TO 600
         IF (X(R).GT.X(R-1)) WNZ = .FALSE.
         IF (W(R).EQ.0.0D+0 .OR. WNZ) GO TO 20
         MDIST = MDIST + 1
         WNZ = .TRUE.
   20 CONTINUE
C
C     CHECK THAT XMIN.LT.XMAX AND THAT XMIN AND XMAX SPAN THE DATA
C     X VALUES.
C
   40 IERROR = 1
      IF (XMIN.GT.X(MFIRST) .OR. XMAX.LT.X(MLAST) .OR. XMIN.GE.XMAX)
     *    GO TO 600
C
C     IF THE NUMBER OF DISTINCT VALUES OF  X(R)  WITH NON-ZERO
C     WEIGHT IS LESS THAN THE NUMBER OF INDEPENDENT COEFFICIENTS
C     IN THE FIT OF MAXIMUM DEGREE  K  THERE IS NO UNIQUE
C     POLYNOMIAL
C     APPROXIMATION OF THAT DEGREE.
C
      L = K - INU
      IERROR = 3
      IF (MDIST.LE.L) GO TO 600
C
C     CHECK THAT  NROWS  HAS BEEN SET SUFFICIENTLY LARGE.
C
      IERROR = 5
      IF (KALL.EQ.1 .AND. NROWS.LT.KPLUS1) GO TO 600
      IF (INUP1.EQ.1) GO TO 80
C
C     NORMALIZE THE FORCING FACTOR SO THAT ITS LEADING COEFFICIENT
C     IS UNITY, CHECKING THAT THIS COEFFICIENT WAS NOT ZERO.
C
      IERROR = 4
      DI = NU(INUP1)
      IF (DI.EQ.0.0D0) GO TO 600
      DO 60 I = 1, INUP1
         WORK2(1,I) = NU(I)/DI
         WORK2(2,I) = 0.0D0
   60 CONTINUE
C
   80 IERROR = 0
      X1 = XMIN
      XM = XMAX
      D = XM - X1
C
C     THE INITIAL VALUES OF EPS(L,R) (L = 1,2,....NDV AND R =
C     MFIRST, MFIRST+1,....MLAST) OF THE WEIGHTED RESIDUALS AND
C     THE VALUES WORK1(1,R)(R=1,2...M) OF THE NORMALIZED
C     INDEPENDENT VARIABLE ARE COMPUTED. N.B. WORK1(1,R) IS
C     COMPUTED FROM THE EXPRESSION BELOW RATHER THAN THE MORE
C     NATURAL FORM   (2.0*X(R) - X1 - XM)/D
C     SINCE THE FORMER GUARANTEES THE COMPUTED VALUE TO DIFFER FROM
C     THE TRUE VALUE BY AT MOST  4.0*MACHINE ACCURACY,  WHEREAS THE
C     LATTER HAS NO SUCH GUARANTEE.
C
C     MDIST IS NOW USED TO RECORD THE TOTAL NUMBER OF DATA POINTS
C     WITH NON-ZERO WEIGHT.
C
      MDIST = 0
      DO 120 R = MFIRST, MLAST
         WR = W(R)
         IF (WR.NE.0.0D0) MDIST = MDIST + 1
         MR = R - MFIRST + 1
         DO 100 L = 1, NDV
            EPS(L,R) = WR*Y(L,R)
  100    CONTINUE
         WORK1(1,MR) = ((X(R)-X1)-(XM-X(R)))/D
  120 CONTINUE
      IM1 = INU*KALL + 1
      BETAI = 0.0D0
      DO 160 JPLUS1 = 1, KPLUS1
         SERR(JPLUS1) = 0.0D0
         DO 140 L = 1, NDV
            A(L,IM1,JPLUS1) = 0.0D0
  140    CONTINUE
  160 CONTINUE
      DO 560 IPLUS1 = INUP1, KPLUS1
C
C        SET STARTING VALUES FOR DEGREE  I.
C
         II = (IPLUS1-1)*KALL + 1
         IPLUS2 = IPLUS1 + 1
         IF (IPLUS1.EQ.KPLUS1) GO TO 240
         IF (KALL.EQ.0) GO TO 220
         DO 200 JPLUS1 = IPLUS2, KPLUS1
            DO 180 L = 1, NDV
               A(L,II,JPLUS1) = 0.0D0
  180       CONTINUE
  200    CONTINUE
  220    WORK2(1,IPLUS2) = 0.0D0
         WORK2(2,IPLUS2) = 0.0D0
  240    ALPIP1 = 0.0D0
         DI = 0.0D0
         DO 260 L = 1, NDV
            CI(L) = 0.0D0
  260    CONTINUE
         WORK2(1,IPLUS1) = 1.0D0
         IF (KPLUS1.GT.1) WORK2(2,1) = WORK2(1,2)
         DO 440 R = MFIRST, MLAST
            IF (W(R).EQ.0.0D0) GO TO 440
            MR = R - MFIRST + 1
            XCAPR = WORK1(1,MR)
C
C           THE WEIGHTED VALUE WORK1(2, R)  OF THE ORTHOGONAL POLYNOMIAL
C           OF DEGREE I AT X = X(R) IS COMPUTED BY RECURRENCE FROM ITS
C           CHEBYSHEV-SERIES REPRESENTATION.
C
            IF (IPLUS1.GT.1) GO TO 280
            WRPR = W(R)*0.5D0*WORK2(1,1)
            WORK1(2,MR) = WRPR
            GO TO 400
  280       J = IPLUS2
            IF (XCAPR.GT.0.5D0) GO TO 360
            IF (XCAPR.GE.-0.5D0) GO TO 320
C
C           GENTLEMAN*S MODIFIED RECURRENCE.
C
            FACTOR = 2.0D0*(1.0D0+XCAPR)
            DJ = 0.0D0
            BJ = 0.0D0
            DO 300 JREV = 2, IPLUS1
               J = J - 1
               DJ = WORK2(1,J) - DJ + FACTOR*BJ
               BJ = DJ - BJ
  300       CONTINUE
            WRPR = W(R)*(0.5D0*WORK2(1,1)-DJ+0.5D0*FACTOR*BJ)
            WORK1(2,MR) = WRPR
            GO TO 400
C
C           CLENSHAW*S ORIGINAL RECURRENCE.
C
  320       FACTOR = 2.0D0*XCAPR
            BJP1 = 0.0D0
            BJ = 0.0D0
            DO 340 JREV = 2, IPLUS1
               J = J - 1
               BJP2 = BJP1
               BJP1 = BJ
               BJ = WORK2(1,J) - BJP2 + FACTOR*BJP1
  340       CONTINUE
            WRPR = W(R)*(0.5D0*WORK2(1,1)-BJP1+0.5D0*FACTOR*BJ)
            WORK1(2,MR) = WRPR
            GO TO 400
C
C           REINSCH*S MODIFIED RECURRENCE.
C
  360       FACTOR = 2.0D0*(1.0D0-XCAPR)
            DJ = 0.0D0
            BJ = 0.0D0
            DO 380 JREV = 2, IPLUS1
               J = J - 1
               DJ = WORK2(1,J) + DJ - FACTOR*BJ
               BJ = BJ + DJ
  380       CONTINUE
            WRPR = W(R)*(0.5D0*WORK2(1,1)+DJ-0.5D0*FACTOR*BJ)
            WORK1(2,MR) = WRPR
C
C           THE COEFFICIENTS CI(L) OF THE I TH ORTHOGONAL POLYNOMIAL
C           L=1,2....NDV AND THE COEFFICIENTS ALPIP1 AND BETAI IN THE
C           THREE-TERM RECURRENCE RELATION FOR THE ORTHOGONAL
C           POLYNOMIALS ARE COMPUTED.
C
  400       WRPRSQ = WRPR**2
            DI = DI + WRPRSQ
            DO 420 L = 1, NDV
               CI(L) = CI(L) + WRPR*EPS(L,R)
  420       CONTINUE
            ALPIP1 = ALPIP1 + WRPRSQ*XCAPR
  440    CONTINUE
         DO 460 L = 1, NDV
            CI(L) = CI(L)/DI
  460    CONTINUE
         IF (IPLUS1.NE.INUP1) BETAI = DI/DIM1
         ALPIP1 = 2.0D0*ALPIP1/DI
C
C        THE WEIGHTED RESIDUALS EPS(L,R)(L=1,2....NDV AND R=MFIRST,
C        MFIRST+1....MLAST) FOR DEGREE I ARE COMPUTED, TOGETHER
C        WITH THEIR SUM OF SQUARES, SIGMAI
C
         DF = MDIST - (IPLUS1-INU)
         DO 500 L = 1, NDV
            CIL = CI(L)
            SIGMAI = 0.0D0
            DO 480 R = MFIRST, MLAST
               IF (W(R).EQ.0.0D0) GO TO 480
               MR = R - MFIRST + 1
               EPSLR = EPS(L,R) - CIL*WORK1(2,MR)
               EPS(L,R) = EPSLR
               SIGMAI = SIGMAI + EPSLR**2
  480       CONTINUE
C
C           THE ROOT MEAN SQUARE RESIDUAL  S(L, I + 1)  FOR DEGREE  I
C           IS THEORETICALLY UNDEFINED IF  M = I + 1 - INU  (THE
C           CONDITION FOR THE POLYNOMIAL TO PASS EXACTLY THROUGH THE
C           DATA POINTS). SHOULD THIS CASE ARISE THE R.M.S. RESIDUAL
C           IS SET TO ZERO.
C
            IF (DF.LE.0.0D0) S(L,IPLUS1) = 0.0D0
            IF (DF.GT.0.0D0) S(L,IPLUS1) = SQRT(SIGMAI/DF)
  500    CONTINUE
C
C        THE CHEBYSHEV COEFFICIENTS A(L, I+1, 1), A(L, I+1, 2)....
C        A(L, I+1, I+1) IN THE POLYNOMIAL APPROXIMATION OF DEGREE I
C        TO EACH SET OF VALUES OF THE INDEPENDENT VARIABLE
C        (L=1,2,...,NDV) TOGETHER WITH THE COEFFICIENTS
C        WORK2(1, 1), WORK2(1, 2), ..., WORK2(1, I + 1),   IN THE
C        CHEBYSHEV-SERIES REPRESENTATION OF THE  (I + 1) TH
C        ORTHOGONAL POLYNOMIAL ARE COMPUTED.
C
         DO 540 JPLUS1 = 1, IPLUS1
            JPLUS2 = JPLUS1 + 1
            PIJ = WORK2(1,JPLUS1)
            SERR(JPLUS1) = SERR(JPLUS1) + PIJ**2/DI
            DO 520 L = 1, NDV
               A(L,II,JPLUS1) = A(L,IM1,JPLUS1) + CI(L)*PIJ
  520       CONTINUE
            IF (JPLUS1.EQ.KPLUS1) GO TO 560
            WORK2(1,JPLUS1) = WORK2(1,JPLUS2) + WORK2(2,JPLUS1) -
     *                        ALPIP1*PIJ - BETAI*WORK2(2,JPLUS2)
            WORK2(2,JPLUS2) = PIJ
  540    CONTINUE
         DIM1 = DI
         IM1 = II
  560 CONTINUE
      DO 580 IPLUS1 = 1, KPLUS1
         SERR(IPLUS1) = 1.0D0/SQRT(SERR(IPLUS1))
  580 CONTINUE
  600 IFAIL = IERROR
      RETURN
      END

      SUBROUTINE E02AEF(NPLUS1,A,XCAP,P,IFAIL)
C     NAG LIBRARY SUBROUTINE  E02AEF
C
C     E02AEF  EVALUATES A POLYNOMIAL FROM ITS CHEBYSHEV-
C     SERIES REPRESENTATION.
C
C     CLENSHAW METHOD WITH MODIFICATIONS DUE TO REINSCH
C     AND GENTLEMAN.
C
C     USES NAG LIBRARY ROUTINES  P01ABF  AND  X02AJF.
C     USES INTRINSIC FUNCTION  ABS.
C
C     STARTED - 1973.
C     COMPLETED - 1976.
C     AUTHOR - MGC AND JGH.
C
C     NAG COPYRIGHT 1975
C     MARK 5 RELEASE
C     MARK 7 REVISED IER-140 (DEC 1978)
C     MARK 9 REVISED. IER-352 (SEP 1981)
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C     MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (APR 1988).
C
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='E02AEF')
C     .. Scalar Arguments ..
      DOUBLE PRECISION  P, XCAP
      INTEGER           IFAIL, NPLUS1
C     .. Array Arguments ..
      DOUBLE PRECISION  A(NPLUS1)
C     .. Local Scalars ..
      DOUBLE PRECISION  BK, BKP1, BKP2, DK, ETA, FACTOR
      INTEGER           IERROR, K, KREV, N, NPLUS2
C     .. Local Arrays ..
      CHARACTER*1       P01REC(1)
C     .. External Functions ..
      DOUBLE PRECISION  X02AJF
      INTEGER           P01ABF
      EXTERNAL          X02AJF, P01ABF
C     .. Intrinsic Functions ..
      INTRINSIC         ABS
C     .. Executable Statements ..
      IERROR = 0
      ETA = X02AJF()
C     INSERT CALL TO X02AJF
C
C     ETA  IS THE SMALLEST POSITIVE NUMBER SUCH THAT
C     THE COMPUTED VALUE OF  1.0 + ETA  EXCEEDS UNITY.
C
      IF (NPLUS1.GE.1) GO TO 20
      IERROR = 2
      GO TO 180
   20 IF (ABS(XCAP).LE.1.0D0+4.0D0*ETA) GO TO 40
      IERROR = 1
      P = 0.0D0
      GO TO 180
   40 IF (NPLUS1.GT.1) GO TO 60
      P = 0.5D0*A(1)
      GO TO 180
   60 N = NPLUS1 - 1
      NPLUS2 = N + 2
      K = NPLUS2
      IF (XCAP.GT.0.5D0) GO TO 140
      IF (XCAP.GE.-0.5D0) GO TO 100
C
C     GENTLEMAN*S MODIFIED RECURRENCE.
C
      FACTOR = 2.0D0*(1.0D0+XCAP)
      DK = 0.0D0
      BK = 0.0D0
      DO 80 KREV = 1, N
         K = K - 1
         DK = A(K) - DK + FACTOR*BK
         BK = DK - BK
   80 CONTINUE
      P = 0.5D0*A(1) - DK + 0.5D0*FACTOR*BK
      GO TO 180
C
C     CLENSHAW*S ORIGINAL RECURRENCE.
C
  100 FACTOR = 2.0D0*XCAP
      BKP1 = 0.0D0
      BK = 0.0D0
      DO 120 KREV = 1, N
         K = K - 1
         BKP2 = BKP1
         BKP1 = BK
         BK = A(K) - BKP2 + FACTOR*BKP1
  120 CONTINUE
      P = 0.5D0*A(1) - BKP1 + 0.5D0*FACTOR*BK
      GO TO 180
C
C     REINSCH*S MODIFIED RECURRENCE.
C
  140 FACTOR = 2.0D0*(1.0D0-XCAP)
      DK = 0.0D0
      BK = 0.0D0
      DO 160 KREV = 1, N
         K = K - 1
         DK = A(K) + DK - FACTOR*BK
         BK = BK + DK
  160 CONTINUE
      P = 0.5D0*A(1) + DK - 0.5D0*FACTOR*BK
  180 IF (IERROR) 200, 220, 200
  200 IFAIL = P01ABF(IFAIL,IERROR,SRNAME,0,P01REC)
      RETURN
  220 IFAIL = 0
      RETURN
      END

      SUBROUTINE E02AFF(NPLUS1,F,A,IFAIL)
C     NAG LIBRARY SUBROUTINE  E02AFF
C
C     E02AFF  COMPUTES THE COEFFICIENTS OF A POLYNOMIAL,
C     IN ITS CHEBYSHEV-SERIES FORM, WHICH INTERPOLATES
C     (PASSES EXACTLY THROUGH) DATA AT A SPECIAL SET OF
C     POINTS.  LEAST-SQUARES POLYNOMIAL APPROXIMATIONS
C     CAN ALSO BE OBTAINED.
C
C     CLENSHAW METHOD WITH MODIFICATIONS DUE TO REINSCH
C     AND GENTLEMAN.
C
C     USES NAG LIBRARY ROUTINES  P01AAF  AND  X01AAF.
C     USES BASIC EXTERNAL FUNCTION  SIN.
C
C     STARTED - 1973.
C     COMPLETED - 1976.
C     AUTHOR - MGC AND JGH.
C
C     NAG COPYRIGHT 1975
C     MARK 5 RELEASE
C     MARK 5B REVISED  IER-73
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='E02AFF')
C     .. Scalar Arguments ..
      INTEGER           IFAIL, NPLUS1
C     .. Array Arguments ..
      DOUBLE PRECISION  A(NPLUS1), F(NPLUS1)
C     .. Local Scalars ..
      DOUBLE PRECISION  BK, BKP1, BKP2, DK, F0, FACTOR, FLI, FLN,
     *                  HALFFN, PI, PIBY2N
      INTEGER           I, IERROR, IPLUS1, J, K, KREV, N, N2, NLESS1
C     .. Local Arrays ..
      CHARACTER*1       P01REC(1)
C     .. External Functions ..
      DOUBLE PRECISION  X01AAF
      INTEGER           P01ABF
      EXTERNAL          X01AAF, P01ABF
C     .. Intrinsic Functions ..
      INTRINSIC         SIN
C     .. Executable Statements ..
      IERROR = 0
      IF (NPLUS1.GT.2) GO TO 40
      IF (NPLUS1.EQ.2) GO TO 20
      IERROR = 1
      GO TO 180
   20 A(1) = F(1) + F(2)
      A(2) = 0.5D0*(F(1)-F(2))
      GO TO 180
C
C     SET THE VALUE OF  PI. INSERT CALL TO X01AAF
C
   40 PI = X01AAF(PI)
      N = NPLUS1 - 1
      FLN = N
      N2 = 2*N
      NLESS1 = N - 1
      PIBY2N = 0.5D0*PI/FLN
      F0 = F(1)
      HALFFN = 0.5D0*F(NPLUS1)
      DO 160 IPLUS1 = 1, NPLUS1
         I = IPLUS1 - 1
         K = NPLUS1
         J = 3*I
         IF (J.GT.N2) GO TO 120
         IF (J.GE.N) GO TO 80
C
C        REINSCH*S MODIFIED RECURRENCE.
C
         FLI = I
         FACTOR = 4.0D0*(SIN(PIBY2N*FLI))**2
         DK = HALFFN
         BK = HALFFN
         DO 60 KREV = 1, NLESS1
            K = K - 1
            DK = F(K) + DK - FACTOR*BK
            BK = BK + DK
   60    CONTINUE
         A(IPLUS1) = (F0+2.0D0*DK-FACTOR*BK)/FLN
         GO TO 160
C
C        CLENSHAW*S ORIGINAL RECURRENCE.
C
   80    FLI = N - 2*I
         FACTOR = 2.0D0*SIN(PIBY2N*FLI)
         BKP1 = 0.0D0
         BK = HALFFN
         DO 100 KREV = 1, NLESS1
            K = K - 1
            BKP2 = BKP1
            BKP1 = BK
            BK = F(K) - BKP2 + FACTOR*BKP1
  100    CONTINUE
         A(IPLUS1) = (F0-2.0D0*BKP1+FACTOR*BK)/FLN
         GO TO 160
C
C        GENTLEMAN*S MODIFIED RECURRENCE.
C
  120    FLI = N - I
         FACTOR = 4.0D0*(SIN(PIBY2N*FLI))**2
         DK = HALFFN
         BK = HALFFN
         DO 140 KREV = 1, NLESS1
            K = K - 1
            DK = F(K) - DK + FACTOR*BK
            BK = DK - BK
  140    CONTINUE
         A(IPLUS1) = (F0-2.0D0*DK+FACTOR*BK)/FLN
  160 CONTINUE
      A(NPLUS1) = 0.5D0*A(NPLUS1)
  180 IF (IERROR) 200, 220, 200
  200 IFAIL = P01ABF(IFAIL,IERROR,SRNAME,0,P01REC)
      RETURN
  220 IFAIL = 0
      RETURN
      END

      SUBROUTINE E02AGF(M,KPLUS1,NROWS,XMIN,XMAX,X,Y,W,MF,XF,YF,LYF,IP,
     *                  A,S,NP1,WRK,LWRK,IWRK,LIWRK,IFAIL)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 9 REVISED. IER-314 (SEP 1981).
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C     MARK 14C REVISED. IER-878 (NOV 1990).
C
C     **************************************************************
C
C     NPL ALGORITHMS LIBRARY ROUTINE CONFIT
C
C     CREATED  20/8/1979      UPDATED  16/5/80     RELEASE NO. 00/05
C
C     AUTHOR... GERALD T ANTHONY.
C     NATIONAL PHYSICAL LABORATORY,
C     TEDDINGTON, MIDDLESEX TW11 0LW,
C     ENGLAND.
C
C     **************************************************************
C
C     E02AGF CALLS E01AEW TO DETERMINE A POLYNOMIAL MU(X) WHICH
C     INTERPOLATES THE GIVEN CONSTRAINTS AND A POLYNOMIAL NU(X)
C     WHICH HAS VALUE ZERO WHERE A CONSTRAINED VALUE IS SPECIFIED.
C     IT THEN CALLS E02ADZ TO FIT Y-MU(X) AS A POLYNOMIAL IN X
C     WITH FACTOR NU(X). FINALLY THE COEFFICIENTS OF MU ARE ADDED
C     TO THESE FITS TO GIVE THE COEFFICIENTS OF THE CONSTRAINED
C     FITS TO Y. ALL POLYNOMIALS ARE EXPRESSED IN CHEBYSHEV
C     SERIES FORM
C
C     INPUT PARAMETERS
C        M        THE NUMBER OF DATA POINTS TO BE FITTED
C        KPLUS1   FITS WITH UP TO KPLUS1 COEFFICIENTS ARE REQUIRED
C        NROWS    FIRST DIMENSION OF ARRAY A WHERE COEFFICIENTS ARE
C                 TO BE STORED
C        XMIN,    END POINTS OF THE RANGE OF THE
C        XMAX     INDEPENDENT VARIABLE
C        X, Y, W  ARRAYS OF DATA VALUES OF THE INDEPENDENT VARIABLE,
C                 DEPENDENT VARIABLE AND WEIGHT, RESPECTIVELY
C        MF       NUMBER OF X VALUES AT WHICH A CONSTRAINT
C                 IS SPECIFIED
C        XF       ARRAY OF VALUES OF THE INDEPENDENT
C                 VARIABLE AT WHICH CONSTRAINTS ARE
C                 SPECIFIED
C        YF       ARRAY OF SPECIFIED VALUES AND DERIVATIVES OF THE
C                 DEPENDENT VARIABLE IN THE ORDER
C                    Y1, Y1 DERIV, Y1 2ND DERIV,...., Y2,....
C        LYF      DIMENSION OF ARRAY YF
C        IP       INTEGER ARRAY OF DEGREES OF DERIVATIVES
C                 SPECIFIED AT EACH POINT XF
C
C     OUTPUT PARAMETERS
C        A        ON EXIT, 2 PARAMETER ARRAY CONTAINING THE
C                 COEFFICIENTS OF THE CHEBYSHEV SERIES
C                 REPRESENTATION OF THE FITS, A(I+1, J+1)
C                 CONTAINS THE COEFFICIENT OF TJ IN THE FIT
C                 OF DEGREE I, I = N,N+1,...,K,  J =
C                 0,1,...,I  WHERE N = NP1 - 1
C        S        ON EXIT, ARRAY CONTAINING THE R.M.S. RESIDUAL FOR
C                 EACH DEGREE OF FIT FROM N TO K
C        NP1      ON EXIT, CONTAINS N + 1, WHERE N IS THE
C                 TOTAL NUMBER OF INTERPOLATION CONDITIONS
C
C        IFAIL    FAILURE INDICATOR
C                  0 - SUCCESSFUL TERMINATION
C                  1 - AT LEAST ONE OF THE FOLLOWING CONDITIONS
C                      HAS BEEN VIOLATED
C                      LYF    AT LEAST N
C                      LWRK   AT LEAST 2*N + 2 + THE LARGER OF
C                             4*M + 3*KPLUS1 AND 8*NP1 +
C                             5*IMAX + MF - 3 WHERE IMAX =
C                             1 + MAX(IP(I))
C                      LIWRK  AT LEAST 2*MF + 2
C                      KPLUS1 AT LEAST NP1
C                      M      AT LEAST 1
C                      NROWS  AT LEAST KPLUS1
C                      MF     AT LEAST 1
C                  2 - FOR SOME I, IP(I) IS LESS THAN 0
C                  3 - XMIN IS NOT STRICTLY LESS THAN XMAX
C                      OR FOR SOME I, XF(I) IS NOT IN RANGE
C                      XMIN TO XMAX OR THE XF(I) ARE NOT
C                      DISTINCT
C                  4 - FOR SOME I, X(I) IS NOT IN RANGE XMIN TO XMAX
C                  5 - THE X(I) ARE NOT NON-DECREASING
C                  6 - THE NUMBER OF DISTINCT VALUES OF X(I) WITH
C                      NON-ZERO WEIGHT IS LESS THAN KPLUS1 - NP1
C                  7 - E01AEW HAS FAILED TO CONVERGE, IE
C                      THE CONSTRAINT CANNOT BE SATISFIED
C                      WITH SUFFICIENT ACCURACY
C
C     WORKSPACE PARAMETERS
C        WRK      REAL WORKSPACE ARRAY
C        LWRK     DIMENSION OF WRK.   LWRK MUST BE AT LEAST
C                 2*N + 2 + THE LARGER OF
C                 4*M + 3*KPLUS1 AND 8*NP1 + 5*IMAX + MF - 3
C                 WHERE IMAX = 1 + MAX(IP(I))
C        IWRK     INTEGER WORKSPACE ARRAY
C        LIWRK    DIMENSION OF IWRK.   LIWRK MUST BE AT LEAST
C                 2*MF + 2
C
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='E02AGF')
C     .. Scalar Arguments ..
      DOUBLE PRECISION  XMAX, XMIN
      INTEGER           IFAIL, KPLUS1, LIWRK, LWRK, LYF, M, MF, NP1,
     *                  NROWS
C     .. Array Arguments ..
      DOUBLE PRECISION  A(NROWS,KPLUS1), S(KPLUS1), W(M), WRK(LWRK),
     *                  X(M), XF(MF), Y(M), YF(LYF)
      INTEGER           IP(MF), IWRK(LIWRK)
C     .. Local Scalars ..
      DOUBLE PRECISION  AMUJ, ONE, XI, XMU, ZERO
      INTEGER           I, IERROR, IM1, IMAX, IPI, IYMUX, J, LW, MDIST,
     *                  N, NANU, NEPS, NSER, NWRK, NWRK1, NWRK2
C     .. Local Arrays ..
      CHARACTER*1       P01REC(1)
C     .. External Functions ..
      INTEGER           P01ABF
      EXTERNAL          P01ABF
C     .. External Subroutines ..
      EXTERNAL          E01AEW, E02ADZ, E02AKF
C     .. Data statements ..
      DATA              ONE, ZERO/1.0D+0, 0.0D+0/
C     .. Executable Statements ..
      IERROR = 1
      IF (MF.LT.1) GO TO 280
      IERROR = 2
      IMAX = 0
      NP1 = 1
      DO 20 I = 1, MF
         IPI = IP(I) + 1
         IF (IPI.LT.1) GO TO 280
         IF (IPI.GT.IMAX) IMAX = IPI
         NP1 = NP1 + IPI
   20 CONTINUE
      N = NP1 - 1
      IERROR = 1
      IF (LYF.LT.N .OR. LIWRK.LT.2*MF+2) GO TO 280
      I = 4*M + 3*KPLUS1
      LW = 8*NP1 + 5*IMAX + MF - 3
      IF (LW.LT.I) LW = I
      LW = LW + 2*NP1
      IF (LW.GT.LWRK .OR. NP1.GT.KPLUS1 .OR. M.LT.1) GO TO 280
      NANU = NP1 + 1
      NWRK = NANU + NP1
      IERROR = 3
      IF (XMAX.LE.XMIN) GO TO 280
      IF (XF(1).GT.XMAX .OR. XF(1).LT.XMIN) GO TO 280
      IF (MF.EQ.1) GO TO 80
      DO 60 I = 2, MF
         XI = XF(I)
         IF (XI.GT.XMAX .OR. XI.LT.XMIN) GO TO 280
         IM1 = I - 1
         DO 40 J = 1, IM1
            IF (XI.EQ.XF(J)) GO TO 280
   40    CONTINUE
   60 CONTINUE
   80 IERROR = 6
      XMU = XMIN
      IF (X(1).EQ.XMIN) XMU = XMAX
      MDIST = 0
      DO 140 I = 1, M
         XI = X(I)
         IF (XI.EQ.XMU .OR. W(I).EQ.ZERO) GO TO 140
         DO 100 J = 1, MF
            IF (XI.EQ.XF(J)) GO TO 120
  100    CONTINUE
         MDIST = MDIST + 1
  120    XMU = XI
  140 CONTINUE
      IF (MDIST.LT.KPLUS1-N) GO TO 280
      IWRK(1) = 1
      IERROR = 1
      CALL E01AEW(MF,XMIN,XMAX,XF,YF,IP,N,NP1,5,20,WRK,WRK(NANU),
     *            WRK(NWRK),LWRK-NWRK+1,IWRK,LIWRK,IERROR)
      IF (IERROR.EQ.0) GO TO 160
      IERROR = 7
      GO TO 280
  160 DO 200 I = 1, M
         IERROR = 1
         CALL E02AKF(N,XMIN,XMAX,WRK,1,NP1,X(I),XMU,IERROR)
         IF (IERROR.EQ.0) GO TO 180
         IERROR = 4
         GO TO 280
  180    IYMUX = NWRK + I - 1
C
C        STORE Y - MU(X) AT ITH DATA POINT
C
         WRK(IYMUX) = Y(I) - XMU
  200 CONTINUE
      NWRK1 = NWRK + M
      NWRK2 = NWRK1 + 2*M
      NSER = NWRK2 + 2*KPLUS1
      NEPS = NSER + KPLUS1
      IERROR = 1
      CALL E02ADZ(1,M,M,KPLUS1,NROWS,1,1,X,WRK(NWRK),W,XMIN,XMAX,NP1,
     *            WRK(NANU),WRK(NWRK1),WRK(NWRK2),A,S,WRK(NSER),
     *            WRK(NEPS),IERROR)
      IF (IERROR.EQ.0) GO TO 220
      IERROR = IERROR + 3
      IF (IERROR.EQ.8) IERROR = 1
      GO TO 280
  220 DO 260 J = 1, N
         AMUJ = WRK(J)
         DO 240 I = NP1, KPLUS1
            A(I,J) = A(I,J) + AMUJ
  240    CONTINUE
  260 CONTINUE
  280 IFAIL = P01ABF(IFAIL,IERROR,SRNAME,0,P01REC)
      RETURN
C     END E02AGF
      END
      SUBROUTINE E02AHF(NP1,XMIN,XMAX,A,IA1,LA,PATM1,ADIF,IADIF1,LADIF,
     *                  IFAIL)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     NPL DATA FITTING LIBRARY ROUTINE CHBDIF
C
C     CREATED 1/5/79    UPDATED 23/1/80     RELEASE NO. 00/03
C
C     AUTHORS.. GERALD T ANTHONY, MAURICE G COX, J GEOFFREY HAYES.
C     NATIONAL PHYSICAL LABORATORY
C     TEDDINGTON, MIDDLESEX, ENGLAND.
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     INPUT PARAMETERS
C        NP1    = N+1 WHERE N IS DEGREE OF GIVEN POLYNOMIAL
C        XMIN   LOWER LIMIT OF RANGE OF X
C        XMAX   UPPER LIMIT OF RANGE OF X
C        A      COEFFICIENTS A0, A1,...AN OF THE GIVEN POLYNOMIAL
C        IA1       ARE STORED IN ARRAY A IN POSITIONS 1, 1+IA1,...
C                  1+N*IA1, RESPECTIVELY
C        LA     THE DECLARED DIMENSION OF ARRAY A
C
C     OUTPUT PARAMETERS
C        PATM1  THE VALUE OF THE GIVEN POLYNOMIAL AT XMIN
C        ADIF   THE COEFFICIENTS OF THE DERIVATIVE POLYNOMIAL
C        IADIF1    ARE RETURNED IN ARRAY ADIF IN POSITIONS
C                  1, 1+IADIF1,...1+(N-1)*IADIF1
C        LADIF  THE DECLARED DIMENSION OF ARRAY ADIF
C        IFAIL  ERROR INDICATOR
C
C     DIFFERENTIATE THE SERIES WITH COEFFICIENTS A OF DEGREE N
C     (I.E. NP1 COEFFICIENTS) TO OBTAIN THE SERIES WITH COEFFICIENTS
C     ADIF OF DEGREE N-1. ALSO SET NEXT HIGHER COEFFICIENT TO ZERO.
C
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='E02AHF')
C     .. Scalar Arguments ..
      DOUBLE PRECISION  PATM1, XMAX, XMIN
      INTEGER           IA1, IADIF1, IFAIL, LA, LADIF, NP1
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LA), ADIF(LADIF)
C     .. Local Scalars ..
      INTEGER           IERROR, N
C     .. Local Arrays ..
      CHARACTER*1       P01REC(1)
C     .. External Functions ..
      INTEGER           P01ABF
      EXTERNAL          P01ABF
C     .. External Subroutines ..
      EXTERNAL          E02AHZ
C     .. Executable Statements ..
      N = NP1 - 1
      IERROR = 1
      IF ( .NOT. (NP1.GE.1 .AND. XMAX.GT.XMIN .AND. IA1.GE.1 .AND.
     *    LA.GT.N*IA1 .AND. IADIF1.GE.1 .AND. LADIF.GT.N*IADIF1))
     *    GO TO 20
      IERROR = 0
      CALL E02AHZ(NP1,XMIN,XMAX,A,IA1,LA,PATM1,ADIF,IADIF1,LADIF)
   20 IFAIL = P01ABF(IFAIL,IERROR,SRNAME,0,P01REC)
      RETURN
C     END E02AHF
      END
      SUBROUTINE E02AHZ(NP1,XMIN,XMAX,A,IA1,LA,PATM1,ADIF,IADIF1,LADIF)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     NPL DATA FITTING LIBRARY ROUTINE AUXCDF
C
C     CREATED 1/5/79    UPDATED 23/1/80     RELEASE NO. 00/03
C
C     AUTHORS.. GERALD T ANTHONY, MAURICE G COX, J GEOFFREY HAYES.
C     NATIONAL PHYSICAL LABORATORY
C     TEDDINGTON, MIDDLESEX, ENGLAND.
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     INPUT PARAMETERS
C        NP1    = N+1 WHERE N IS DEGREE OF GIVEN POLYNOMIAL
C        XMIN   LOWER LIMIT OF RANGE OF X
C        XMAX   UPPER LIMIT OF RANGE OF X
C        A      COEFFICIENTS A0, A1,...AN OF THE GIVEN POLYNOMIAL
C        IA1       ARE STORED IN ARRAY A IN POSITIONS 1, 1+IA1,...
C                  1+N*IA1, RESPECTIVELY
C        LA     THE DECLARED DIMENSION OF ARRAY A
C
C     OUTPUT PARAMETERS
C        PATM1  THE VALUE OF THE GIVEN POLYNOMIAL AT XMIN
C        ADIF   THE COEFFICIENTS OF THE DERIVATIVE POLYNOMIAL
C        IADIF1    ARE RETURNED IN ARRAY ADIF IN POSITIONS
C                  1, 1+IADIF1,...1+(N-1)*IADIF1
C        LADIF  THE DECLARED DIMENSION OF ARRAY ADIF
C
C     DIFFERENTIATE THE SERIES WITH COEFFICIENTS A OF DEGREE N
C     (I.E. NP1 COEFFICIENTS) TO OBTAIN THE SERIES WITH COEFFICIENTS
C     ADIF OF DEGREE N-1. ALSO SET NEXT HIGHER COEFFICIENT TO ZERO.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  PATM1, XMAX, XMIN
      INTEGER           IA1, IADIF1, LA, LADIF, NP1
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LA), ADIF(LADIF)
C     .. Local Scalars ..
      DOUBLE PRECISION  PTEMP, R, SCLFTR, TWO, U, V, W
      INTEGER           I, N, NA, NADIF
C     .. Data statements ..
      DATA              TWO/2.0D+0/
C     .. Executable Statements ..
      U = 0.0D+0
      V = U
      SCLFTR = TWO/(XMAX-XMIN)
      N = NP1 - 1
      NADIF = N*IADIF1 + 1
      PTEMP = U
      IF (N.EQ.0) GO TO 40
      NA = N*IA1 + 1
      DO 20 I = 1, N
         R = NP1 - I
         W = U + TWO*R*A(NA)
         PTEMP = A(NA) - PTEMP
C
C        STORE COEFF FORMED PREVIOUS TIME ROUND. FIRST TIME ROUND
C        STORE ZERO AS COEFF OF DEGREE N.
C
         ADIF(NADIF) = SCLFTR*V
         U = V
         V = W
         NA = NA - IA1
         NADIF = NADIF - IADIF1
   20 CONTINUE
   40 ADIF(NADIF) = SCLFTR*V
      PATM1 = A(1)/TWO - PTEMP
      RETURN
C     END OF E02AHZ
      END
      SUBROUTINE E02AJF(NP1,XMIN,XMAX,A,IA1,LA,QATM1,AIN,IAINT1,LAINT,
     *                  IFAIL)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 9 REVISED. IER-315 (SEP 1981).
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     NPL DATA FITTING LIBRARY ROUTINE CHBINT
C
C     CREATED 1/5/79    UPDATED 23/1/80     RELEASE NO. 00/03
C
C     AUTHORS.. GERALD T ANTHONY, MAURICE G COX, J GEOFFREY HAYES.
C     NATIONAL PHYSICAL LABORATORY
C     TEDDINGTON, MIDDLESEX, ENGLAND.
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     INPUT PARAMETERS
C        NP1    = N+1 WHERE N IS DEGREE OF GIVEN POLYNOMIAL
C        XMIN   LOWER LIMIT OF RANGE OF X
C        XMAX   UPPER LIMIT OF RANGE OF X
C        A      COEFFICIENTS A0, A1,...AN OF THE GIVEN POLYNOMIAL
C        IA1       ARE STORED IN ARRAY A IN POSITIONS 1, 1+IA1,...
C                  1+N*IA1, RESPECTIVELY
C        LA     THE DECLARED DIMENSION OF ARRAY A
C        QATM1  THE VALUE OF THE INTEGRATED POLYNOMIAL AT XMIN
C
C     OUTPUT PARAMETERS
C        AIN   THE COEFFICIENTS OF THE INTEGRATED POLYNOMIAL
C        IAINT1    ARE RETURNED IN ARRAY AIN IN POSITIONS
C                  1, 1+IAINT1,...1+NP1*IAINT1
C        LAINT  THE DECLARED DIMENSION OF ARRAY AIN
C        IFAIL  ERROR INDICATOR
C
C     INTEGRATE THE SERIES WITH COEFFICIENTS A OF DEGREE N
C     (I.E. NP1 COEFFICIENTS) TO OBTAIN THE SERIES WITH COEFFICIENTS
C     AIN OF DEGREE N + 1.  THE SUM OF THE INTEGRATED SERIES IS
C     QATM1 AT THE LEFT HAND END OF THE INTERVAL OF DEFINITION, XMIN
C
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='E02AJF')
C     .. Scalar Arguments ..
      DOUBLE PRECISION  QATM1, XMAX, XMIN
      INTEGER           IA1, IAINT1, IFAIL, LA, LAINT, NP1
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LA), AIN(LAINT)
C     .. Local Scalars ..
      INTEGER           IERROR
C     .. Local Arrays ..
      CHARACTER*1       P01REC(1)
C     .. External Functions ..
      INTEGER           P01ABF
      EXTERNAL          P01ABF
C     .. External Subroutines ..
      EXTERNAL          E02AJZ
C     .. Executable Statements ..
      IERROR = 1
      IF ( .NOT. (NP1.GE.1 .AND. XMAX.GT.XMIN .AND. LA.GT.(NP1-1)
     *    *IA1 .AND. IA1.GE.1 .AND. IAINT1.GE.1 .AND. LAINT.GT.NP1*
     *    IAINT1)) GO TO 20
      IERROR = 0
      CALL E02AJZ(NP1,XMIN,XMAX,A,IA1,LA,QATM1,AIN,IAINT1,LAINT)
   20 IFAIL = P01ABF(IFAIL,IERROR,SRNAME,0,P01REC)
      RETURN
C     END OF E02AJF
      END
      SUBROUTINE E02AJZ(NP1,XMIN,XMAX,A,IA1,LA,QATM1,AIN,IAINT1,LAINT)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 9 REVISED. IER-315 (SEP 1981).
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     NPL DATA FITTING LIBRARY ROUTINE AUXCIN
C
C     CREATED 1/5/79    UPDATED 23/1/80     RELEASE NO. 00/03
C
C     AUTHORS.. GERALD T ANTHONY, MAURICE G COX, J GEOFFREY HAYES.
C     NATIONAL PHYSICAL LABORATORY
C     TEDDINGTON, MIDDLESEX, ENGLAND.
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     INPUT PARAMETERS
C        NP1    = N+1 WHERE N IS DEGREE OF GIVEN POLYNOMIAL
C        XMIN   LOWER LIMIT OF RANGE OF X
C        XMAX   UPPER LIMIT OF RANGE OF X
C        A      COEFFICIENTS A0, A1,...AN OF THE GIVEN POLYNOMIAL
C        IA1       ARE STORED IN ARRAY A IN POSITIONS 1, 1+IA1,...
C                  1+N*IA1, RESPECTIVELY
C        LA     THE DECLARED DIMENSION OF ARRAY A
C        QATM1  THE VALUE OF THE INTEGRATED POLYNOMIAL AT XMIN
C
C     OUTPUT PARAMETERS
C        AIN   THE COEFFICIENTS OF THE INTEGRATED POLYNOMIAL
C        IAINT1    ARE RETURNED IN ARRAY AIN IN POSITIONS
C                  1, 1+IAINT1,...1+NP1*IAINT1
C        LAINT  THE DECLARED DIMENSION OF ARRAY AIN
C
C     INTEGRATE THE SERIES WITH COEFFICIENTS A OF DEGREE N
C     (I.E. NP1 COEFFICIENTS) TO OBTAIN THE SERIES WITH COEFFICIENTS
C     AIN OF DEGREE N + 1.  THE SUM OF THE INTEGRATED SERIES IS
C     QATM1 AT THE LEFT HAND END OF THE INTERVAL OF DEFINITION, XMIN
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  QATM1, XMAX, XMIN
      INTEGER           IA1, IAINT1, LA, LAINT, NP1
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LA), AIN(LAINT)
C     .. Local Scalars ..
      DOUBLE PRECISION  AR, ARM, ARP, R, SCLFTR, TWO, ZERO
      INTEGER           I, NA, NAINT
C     .. Data statements ..
      DATA              TWO, ZERO/2.0D+0, 0.0D+0/
C     .. Executable Statements ..
      ARP = ZERO
      AR = ARP
      SCLFTR = (XMAX-XMIN)/TWO
      NA = (NP1-1)*IA1 + 1
      NAINT = NP1*IAINT1 + 1
      DO 20 I = 1, NP1
         R = NP1 - I + 1
         ARM = A(NA)
         AIN(NAINT) = SCLFTR*(ARM-ARP)/(TWO*R)
         ARP = AR
         AR = ARM
         NA = NA - IA1
         NAINT = NAINT - IAINT1
   20 CONTINUE
C
C     FORM CONSTANT COEFF SO THAT THE VALUE OF THE INTEGRATED SERIES
C     IS QATM1 AT THE LOWER END OF THE RANGE OF THE INDEP. VARIABLE.
C
      AR = ZERO
      NAINT = NP1*IAINT1 + 1
      DO 40 I = 1, NP1
         AR = AIN(NAINT) - AR
         NAINT = NAINT - IAINT1
   40 CONTINUE
      AIN(1) = TWO*(QATM1+AR)
      RETURN
C     END OF E02AJZ
      END
      SUBROUTINE E02AKF(NP1,XMIN,XMAX,A,IA1,LA,X,RESULT,IFAIL)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     NPL DATA FITTING LIBRARY ROUTINE TVAL1C
C
C     CREATED 19/3/79    UPDATED 6/7/79    RELEASE NO. 00/04.
C
C     AUTHORS.. GERALD T ANTHONY, MAURICE G COX, BETTY CURTIS
C     AND J GEOFFREY HAYES.
C     NATIONAL PHYSICAL LABORATORY
C     TEDDINGTON, MIDDLESEX, ENGLAND.
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     INPUT PARAMETERS
C        NP1      NP1 = N + 1. N IS THE DEGREE OF THE
C                 CHEBYSHEV SERIES
C        XMIN     MINIMUM VALUE OF X
C        XMAX     MAXIMUM VALUE OF X
C        A        THE ARRAY WHERE THE COEFFICIENTS ARE STORED
C        IA1      THE ADDRESS INCREMENT OF A
C        LA       DIMENSION OF A
C        X        UNNORMALIZED ARGUMENT IN THE RANGE (XMIN, XMAX)
C
C     OUTPUT PARAMETERS
C        RESULT   VALUE OF THE SUMMATION
C        IFAIL    ERROR INDICATOR
C
C     NP1 CHEBYSHEV COEFFICIENTS A0, A1, ..., AN, ARE
C     STORED IN THE ARRAY A IN POSITIONS 1, 1+IA1, 1+2*IA1, ...,
C     1+N*IA1, WHERE N = NP1 - 1.
C     IA1 MUST NOT BE NEGATIVE.
C     LA MUST BE AT LEAST EQUAL TO 1 + N*IA1.
C     THE VALUE OF THE POLYNOMIAL OF DEGREE N
C     A0T0(XCAP)/2 + A1T1(XCAP) + A2T2(XCAP) + + ... + ANTN(XCAP),
C     IS CALCULATED FOR THE ARGUMENT XCAP, WHERE XCAP IS
C     THE NORMALIZED VALUE OF X IN THE RANGE (XMIN, XMAX),
C     STORING IT IN RESULT.
C     UNLESS THE ROUTINE DETECTS AN ERROR, IFAIL CONTAINS
C     ZERO ON EXIT.
C     IFAIL = 1 INDICATES AT LEAST ONE OF THE RESTRICTIONS ON
C        INPUT PARAMETERS IS VIOLATED - IE
C     NP1 .GT. 0
C     IA1 .GE. 0
C     LA .GE. 1 + N * IA1
C     XMIN .LT. XMAX
C     IFAIL = 2 INDICATES THAT
C     X DOES NOT SATISFY THE RESTRICTION XMIN .LE. X .LE. XMAX.
C     THE RECURRENCE RELATION BY CLENSHAW, MODIFIED BY REINSCH
C     AND GENTLEMAN, IS USED.
C
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='E02AKF')
C     .. Scalar Arguments ..
      DOUBLE PRECISION  RESULT, X, XMAX, XMIN
      INTEGER           IA1, IFAIL, LA, NP1
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LA)
C     .. Local Scalars ..
      DOUBLE PRECISION  XCAP
      INTEGER           IERROR
C     .. Local Arrays ..
      CHARACTER*1       P01REC(1)
C     .. External Functions ..
      INTEGER           P01ABF
      EXTERNAL          P01ABF
C     .. External Subroutines ..
      EXTERNAL          E02AKY, E02AKZ
C     .. Executable Statements ..
      IERROR = 1
      IF (NP1.LT.1) GO TO 20
      IF (IA1.LT.1) GO TO 20
      IF (LA.LT.1+(NP1-1)*IA1) GO TO 20
      IF (XMAX.LE.XMIN) GO TO 20
      IERROR = IERROR + 1
      IF ((X.GT.XMAX) .OR. (X.LT.XMIN)) GO TO 20
      IERROR = 0
      CALL E02AKY(XMIN,XMAX,X,XCAP)
      CALL E02AKZ(NP1,A,IA1,LA,XCAP,RESULT)
   20 IFAIL = P01ABF(IFAIL,IERROR,SRNAME,0,P01REC)
      RETURN
C
C     END E02AKF
C
      END
      SUBROUTINE E02AKY(XMIN,XMAX,X,XCAP)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     NPL DATA FITTING LIBRARY ROUTINE NRMLZE
C
C     CREATED 5/5/78    UPDATED 11/12/78    RELEASE NO. 00/02.
C
C     AUTHORS.. GERALD T ANTHONY, MAURICE G COX, BETTY CURTIS
C     AND J GEOFFREY HAYES.
C     NATIONAL PHYSICAL LABORATORY
C     TEDDINGTON, MIDDLESEX, ENGLAND.
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     INPUT PARAMETERS
C        XMIN     MINIMUM VALUE OF X
C        XMAX     MAXIMUM VALUE OF X
C        X        UNNORMALIZED ARGUMENT IN RANGE (XMIN, XMAX)
C
C     OUTPUT PARAMETER
C        XCAP     NORMALIZED VALUE OF X
C
C     A VALUE OF X IS GIVEN, SUCH THAT
C     XMIN .LE. X .LE. XMAX.
C     XCAP IS CALCULATED SO THAT -1 .LE. X .LE. +1.
C
C     THIS FORM FOR XCAP ENSURES THAT THE COMPUTED VALUE HAS A
C     VERY SMALL ABSOLUTE ERROR.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  X, XCAP, XMAX, XMIN
C     .. Executable Statements ..
      XCAP = ((X-XMIN)-(XMAX-X))/(XMAX-XMIN)
      RETURN
C
C     END E02AKY
C
      END
      SUBROUTINE E02AKZ(NP1,A,IA1,LA,XCAP,RESULT)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     NPL DATA FITTING LIBRARY ROUTINE TVAL1
C
C     CREATED 9/5/78    UPDATED 6/4/79    RELEASE NO. 00/07.
C
C     AUTHORS.. GERALD T ANTHONY, MAURICE G COX, BETTY CURTIS
C     AND J GEOFFREY HAYES.
C     NATIONAL PHYSICAL LABORATORY
C     TEDDINGTON, MIDDLESEX, ENGLAND.
C
C     * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C     INPUT PARAMETERS
C        NP1      NP1 = N + 1. N IS THE DEGREE OF THE
C                 CHEBYSHEV SERIES
C        A        THE ARRAY WHERE THE COEFFICIENTS ARE STORED
C        IA1      THE ADDRESS INCREMENT OF A
C        LA       DIMENSION OF A
C        XCAP     NORMALIZED ARGUMENT OF THE POLYNOMIAL
C
C     OUTPUT PARAMETER
C        RESULT   VALUE OF THE SUMMATION
C
C     NP1 CHEBYSHEV COEFFICIENTS A0, A1, ..., AN, ARE
C     STORED IN THE ARRAY A IN POSITIONS 1, 1+IA1, 1+2*IA1, ...,
C     1+N*IA1, WHERE N = NP1 - 1.
C     IA1 MUST NOT BE NEGATIVE.
C     LA MUST BE AT LEAST EQUAL TO 1 + N*IA1.
C     THE ARGUMENT XCAP IS ASSUMED TO LIE IN THE RANGE
C     -1 .LE. XCAP .LE. +1.
C     THE VALUE OF THE POLYNOMIAL OF DEGREE N
C     A0T0(XCAP)/2 + A1T1(XCAP) + A2T2(XCAP) + + ... + ANTN(XCAP),
C     IS CALCULATED FOR THE ARGUMENT XCAP STORING IT IN RESULT.
C     THE RECURRENCE RELATION BY CLENSHAW, MODIFIED BY REINSCH
C     AND GENTLEMAN, IS USED.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  RESULT, XCAP
      INTEGER           IA1, LA, NP1
C     .. Array Arguments ..
      DOUBLE PRECISION  A(LA)
C     .. Local Scalars ..
      DOUBLE PRECISION  AJ, BJ, CJ, FACTOR, HALF, SUM, TWO, ZERO
      INTEGER           J, JREV, N
C     .. Data statements ..
      DATA              ZERO, HALF, TWO/0.0D0, 0.5D0, 2.0D0/
C     .. Executable Statements ..
      IF (NP1.GT.1) GO TO 20
      SUM = HALF*A(1)
      GO TO 140
   20 N = NP1 - 1
      AJ = ZERO
      BJ = ZERO
      J = 1 + NP1*IA1
      IF (XCAP.GT.HALF) GO TO 100
      IF (XCAP.GE.-HALF) GO TO 60
C
C     GENTLEMANS MODIFIED RECURRENCE.
C
      FACTOR = TWO + (XCAP+XCAP)
C
C     BRACKETING NECESSARY SO AS TO AVOID ERRORS
C
      DO 40 JREV = 1, N
         J = J - IA1
         AJ = A(J) - AJ + BJ*FACTOR
         BJ = AJ - BJ
   40 CONTINUE
      SUM = HALF*A(1) - AJ + HALF*FACTOR*BJ
      GO TO 140
C
C     CLENSHAWS ORIGINAL RECURRENCE.
C
   60 FACTOR = XCAP + XCAP
      DO 80 JREV = 1, N
         J = J - IA1
         CJ = BJ
         BJ = AJ
         AJ = A(J) - CJ + BJ*FACTOR
   80 CONTINUE
      SUM = HALF*A(1) - BJ + HALF*FACTOR*AJ
      GO TO 140
C
C     REINSCHS MODIFIED RECURRENCE.
C
  100 FACTOR = TWO - (XCAP+XCAP)
C
C     BRACKETING NECESSARY IN ORDER TO AVOID ERRORS
C
      DO 120 JREV = 1, N
         J = J - IA1
         AJ = A(J) + AJ - BJ*FACTOR
         BJ = AJ + BJ
  120 CONTINUE
      SUM = HALF*A(1) + AJ - HALF*FACTOR*BJ
  140 RESULT = SUM
      RETURN
C
C     END E02AKZ
C
      END
